#include "blurmap.h"
#include "../data/elementInformation.h"
#include "../contour/mcconextractor.h"

using namespace PBVP;

BlurMap::BlurMap()
{
	m_Blobby = -2.3442; // this value is from the following paper
	m_IsoDensity = 1.0f;
	//"Evaluation of Protein Docking Predictions Using Hex 3.1 in CAPRI Rounds 1
	// and 2", D.W. Ritchie (2003) PROTEINS: Struct. Funct. Genet. 52(1), 98-106.
	m_Dim[0] = m_Dim[1] = m_Dim[2] = 96;
	m_MinExt[0] = m_MinExt[1] = m_MinExt[2] = 0;
	m_MaxExt[0] = m_MaxExt[1] = m_MaxExt[2] = 1;
}

Reg3DVolumeData* BlurMap::constructVolume(const QVector<Atom> &atoms)
{
	Reg3DVolumeData *pReg3Data = new Reg3DVolumeData(m_Dim);
	computeMinMaxExt(atoms);
	pReg3Data->setMinMaxExt(m_MinExt, m_MaxExt);
	pReg3Data->setNumOfVariables(1);
	float *data = new float[m_Dim[0]*m_Dim[1]*m_Dim[2]];
	memset(data, 0, sizeof(float)*m_Dim[0]*m_Dim[1]*m_Dim[2]);

	double orig[3], span[3];
	for (int i = 0; i < 3; i++)
	{
		orig[i] = m_MinExt[i];
		span[i] = (m_MaxExt[i] - m_MinExt[i])/(m_Dim[i]-1);
	}
	const double eps = 1.0e-4;
	double maxRadFactor = sqrt(1.0 + log(eps) / m_Blobby);  // why??
	// compute contributions from each atom
	for (int i = 0; i < atoms.size(); i++)
	{
		Atom atom = atoms[i];
		double radius = atom.radius;
		double maxRadius = radius * maxRadFactor;
		int minIndex[3], maxIndex[3];
		for(int j = 0; j < 3; j++) {
			minIndex[j] = (int)((atom.position[j]- maxRadius-m_MinExt[j])/span[j]);
			minIndex[j] = (minIndex[j] < 0)? 0:minIndex[j];
			maxIndex[j] = (int)((atom.position[j]+maxRadius-m_MinExt[j])/span[j]);
			maxIndex[j] = (maxIndex[j] > (m_Dim[j]-1))? (m_Dim[j]-1): maxIndex[j];
		}
		for (int iz = minIndex[2]; iz <= maxIndex[2]; iz++)
			for (int iy = minIndex[1]; iy <= maxIndex[1]; iy++ )
				for(int ix = minIndex[0]; ix <= maxIndex[0]; ix++)
		{
			Vector3 pos(orig[0]+ix*span[0], orig[1]+iy*span[1], orig[2]+iz*span[2]);
			int n = ix + iy*m_Dim[0] + iz*m_Dim[0]*m_Dim[1];
			float val = evalDensity(atom, pos);
			data[n] += val;
		}
	}

	pReg3Data->setDataArray(data);
	return pReg3Data;
}

GeometryData* BlurMap::constructSurface(const QVector<Atom>& atoms)
{
	Reg3DVolumeData* pReg3Data = constructVolume(atoms);
	MCContourExtractor extractor;
	GeometryData *pContour = extractor.extractContour(pReg3Data, m_IsoDensity);
	delete pReg3Data;
	return pContour;
}

float BlurMap::evalDensity(const Atom& atom, const Vector3& pos)
{
	double r_2 = (atom.position.x - pos.x)*(atom.position.x - pos.x) + 
				 (atom.position.y - pos.y)*(atom.position.y - pos.y) +
				 (atom.position.z - pos.z)*(atom.position.z - pos.z);
	double r0_2 = atom.radius * atom.radius;
	double expval = exp(m_Blobby*r_2/r0_2 - m_Blobby);
	return expval;
}

void BlurMap::computeMinMaxExt(const QVector<Atom> &atoms)
{
	m_MinExt[0] = m_MaxExt[0] = atoms[0].position[0];
	m_MinExt[1] = m_MaxExt[1] = atoms[0].position[1];
	m_MinExt[2] = m_MaxExt[2] = atoms[0].position[2];
	for (int i = 1; i < atoms.size(); i++)
	{
		if (m_MinExt[0] > atoms[i].position[0]) m_MinExt[0] = atoms[i].position[0];
		if (m_MinExt[1] > atoms[i].position[1]) m_MinExt[1] = atoms[i].position[1];
		if (m_MinExt[2] > atoms[i].position[2]) m_MinExt[2] = atoms[i].position[2];
		if (m_MaxExt[0] < atoms[i].position[0]) m_MaxExt[0] = atoms[i].position[0];
		if (m_MaxExt[1] < atoms[i].position[1]) m_MaxExt[1] = atoms[i].position[1];
		if (m_MaxExt[2] < atoms[i].position[2]) m_MaxExt[2] = atoms[i].position[2];
	}
	double offset = 4*1.5;
	m_MinExt[0] -= offset; 
	m_MinExt[1] -= offset; 
	m_MinExt[2] -= offset;
	m_MaxExt[0] += offset; 
	m_MaxExt[1] += offset;
	m_MaxExt[2] += offset;
}